
import geopandas as pd
import contextily as ctx # Used for contextual basemaps
import matplotlib.pyplot as plt
from geocube.api.core import make_geocube # Used for rasterizing
import os
import shapely
import imageio # Used for making animated GIFs
import numpy as np
from IPython.display import Image
from osgeo import gdal # Raster operations
import zipfile
import rasterio
import rasterio.merge
import rasterio.plot
import rasterio.mask
plt.rcParams['figure.figsize'] = (20, 20)
os.listdir("input")

/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.0-CAPI-1.16.2). Conversions between both will be slow.
  warnings.warn(

['lds-nz-road-centrelines-topo-150k-FGDB.zip',
 'lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip',
 'statsnzpopulation-by-meshblock-2013-census-FGDB.zip',
 'statsnz2018-census-electoral-population-meshblock-2020-FGDB.zip',
 'statsnzregional-council-2021-clipped-generalised-FGDB.zip',
 'lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip']


%%time
REGC = pd.read_file("input/statsnzregional-council-2021-clipped-generalised-FGDB.zip!regional-council-2021-clipped-generalised.gdb")
AKL = REGC[REGC.REGC2021_V1_00_NAME == "Auckland Region"].copy()
# Filter out islands
AKL["geometry"] = max(AKL.geometry.explode(), key=lambda a: a.area)
# Coordinate reference system (projection)
print(AKL.crs)
# Simplify geometry to speed up clip operations
AKL = AKL.simplify(1000).buffer(1000)
ax = AKL.to_crs(epsg=3857).boundary.plot()
ax.set_title("Auckland Region clip extent")
ctx.add_basemap(ax)

epsg:2193
CPU times: user 1.83 s, sys: 936 ms, total: 2.77 s
Wall time: 14.3 s


%%time
df = pd.read_file("zip://input/lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip!lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gdb")

CPU times: user 1min 43s, sys: 2.2 s, total: 1min 46s
Wall time: 1min 46s


print(df.columns)
print(df.crs)
display(df.sample(5))

Index(['Name_2018', 'Name_2012', 'Name_2008', 'Name_2001', 'Name_1996',
       'Class_2018', 'Class_2012', 'Class_2008', 'Class_2001', 'Class_1996',
       'Wetland_18', 'Wetland_12', 'Wetland_08', 'Wetland_01', 'Wetland_96',
       'Onshore_18', 'Onshore_12', 'Onshore_08', 'Onshore_01', 'Onshore_96',
       'EditAuthor', 'EditDate', 'LCDB_UID', 'geometry'],
      dtype='object')
epsg:2193


%%time
df = pd.clip(df, AKL)

CPU times: user 46.2 s, sys: 10.3 ms, total: 46.3 s
Wall time: 46.2 s


df.sample(5)


df.Name_2018.value_counts()

Exotic Forest                                3981
Indigenous Forest                            3673
Manuka and/or Kanuka                         2282
Broadleaved Indigenous Hardwoods             1788
Built-up Area (settlement)                   1350
High Producing Exotic Grassland              1326
Mangrove                                     1151
Urban Parkland/Open Space                    1099
Estuarine Open Water                          441
Orchard, Vineyard or Other Perennial Crop     436
Short-rotation Cropland                       362
Lake or Pond                                  326
Herbaceous Saline Vegetation                  303
Low Producing Grassland                       291
Gorse and/or Broom                            287
Forest - Harvested                            266
Sand or Gravel                                252
Deciduous Hardwoods                           201
Surface Mine or Dump                          132
Mixed Exotic Shrubland                        120
Herbaceous Freshwater Vegetation              118
Transport Infrastructure                      107
River                                          15
Gravel or Rock                                  9
Flaxland                                        9
Fernland                                        2
Matagouri or Grey Scrub                         1
Name: Name_2018, dtype: int64


def simplify_classes(code):
    if code in [1, 2, 5]:
        return 1, "Urban"
    elif code in [68,69,71]:
        return 2, "Vegetation"
    elif code in [0,20,21,22,45,46]:
        return 3, "Water"
    else:
        return 4, "Other"

summary = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
    print(year)
    class_year = f"Class_{year}"
    df[class_year + "_simplified_code"] = df[class_year].apply(lambda c: simplify_classes(c)[0])
    df[class_year + "_simplified_name"] = df[class_year].apply(lambda c: simplify_classes(c)[1])
    summary.append(df[class_year + "_simplified_name"].value_counts())

1996
2001
2008
2012
2018


pd.GeoDataFrame(summary).plot.area()

<AxesSubplot:>


%%capture
# %%capture suppresses output
ims = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
    ax = df.plot(column=f'Class_{year}_simplified_name', legend=True)
    ax.set_title(year)
    ax.figure.tight_layout()
    canvas = ax.figure.canvas
    canvas.draw() # draw the canvas, cache the renderer
    image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
    image = image.reshape(canvas.get_width_height()[::-1] + (3,))
    ims.append(image)

imageio.mimsave("land_use.gif", ims, fps=1)


with open('land_use.gif','rb') as file:
    display(Image(file.read()))


cols = [f"Class_{year}_simplified_code" for year in years]
cols

['Class_1996_simplified_code',
 'Class_2001_simplified_code',
 'Class_2008_simplified_code',
 'Class_2012_simplified_code',
 'Class_2018_simplified_code']


%%time
geocube = make_geocube(
    vector_data=df,
    output_crs="epsg:2193",
    measurements=cols,
    resolution=(-100, 100),
    fill=0, # NaNs, like offshore areas, will be 0
)
geocube

CPU times: user 20 s, sys: 225 µs, total: 20 s
Wall time: 20 s

<xarray.Dataset>
Dimensions:                     (x: 1001, y: 1320)
Coordinates:
  * y                           (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
  * x                           (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
    spatial_ref                 int64 0
Data variables:
    Class_1996_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2001_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2008_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2012_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2018_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


geocube.Class_2018_simplified_code.plot()

<matplotlib.collections.QuadMesh at 0x7f33c2bcaf40>


for year in years:
    print(year)
    outfile = f"output/land_use_{year}.tif"
    if not os.path.isfile(outfile):
        geocube[f"Class_{year}_simplified_code"].rio.to_raster(outfile, dtype=np.byte) # Use np.byte for smaller output filesize

1996
2001
2008
2012
2018


%%time
pop2013 = pd.read_file("input/statsnzpopulation-by-meshblock-2013-census-FGDB.zip!population-by-meshblock-2013-census.gdb")

CPU times: user 9.05 s, sys: 40.9 ms, total: 9.09 s
Wall time: 9.07 s


%%time
pop2013 = pd.clip(pop2013, AKL)

CPU times: user 18.7 s, sys: 135 µs, total: 18.7 s
Wall time: 18.7 s


display(pop2013.sample(5))


#pop2013.Population_Count_Usual_Resident_2013.replace(0, np.nan, inplace=True)


pop2013.Population_Count_Usual_Resident_2013.plot(kind="hist", bins=200)

<AxesSubplot:ylabel='Frequency'>


%%time
pop2013_cube = make_geocube(
    vector_data=pop2013,
    measurements=["Population_Count_Usual_Resident_2013"],
    like=geocube, # Ensures dimensions match
    fill=0 # NaNs, like offshore areas, will be 0
)
pop2013_cube

CPU times: user 2.75 s, sys: 10.2 ms, total: 2.76 s
Wall time: 2.75 s

<xarray.Dataset>
Dimensions:                               (x: 1001, y: 1320)
Coordinates:
  * y                                     (y) float64 6.002e+06 ... 5.87e+06
  * x                                     (x) float64 1.704e+06 ... 1.804e+06
    spatial_ref                           int64 0
Data variables:
    Population_Count_Usual_Resident_2013  (y, x) float64 0.0 0.0 0.0 ... 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


pop2013_cube.Population_Count_Usual_Resident_2013.plot()
outfile = "output/pop2013.tif"
if not os.path.isfile(outfile):
    # byte max value is 255, and we have larger values than that here. uint16 max value is 65535, which is fine
    pop2013_cube.Population_Count_Usual_Resident_2013.rio.to_raster(outfile, dtype=np.uint16)


%%time
pop2018 = pd.read_file("input/statsnz2018-census-electoral-population-meshblock-2020-FGDB.zip!2018-census-electoral-population-meshblock-2020.gdb")

CPU times: user 10.3 s, sys: 30.4 ms, total: 10.3 s
Wall time: 10.3 s


%%time
pop2018 = pd.clip(pop2018, AKL)

CPU times: user 21.3 s, sys: 0 ns, total: 21.3 s
Wall time: 21.3 s


display(pop2018.sample(5))


pop2018.General_Electoral_Population.replace(-999, 0, inplace=True)


pop2018.General_Electoral_Population.plot(kind="hist", bins=200)

<AxesSubplot:ylabel='Frequency'>


%%time
pop2018_cube = make_geocube(
    vector_data=pop2018,
    measurements=["General_Electoral_Population"],
    like=geocube, # Ensures dimensions match
    fill=0
)
pop2018_cube

CPU times: user 2.47 s, sys: 9.71 ms, total: 2.48 s
Wall time: 2.47 s

<xarray.Dataset>
Dimensions:                       (x: 1001, y: 1320)
Coordinates:
  * y                             (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
  * x                             (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
    spatial_ref                   int64 0
Data variables:
    General_Electoral_Population  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


pop2018_cube.General_Electoral_Population.plot()
outfile = "output/pop2018.tif"
if not os.path.isfile(outfile):
    pop2018_cube.General_Electoral_Population.rio.to_raster(outfile, dtype=np.uint16)


%%time
roads = pd.read_file("input/lds-nz-road-centrelines-topo-150k-FGDB.zip!nz-road-centrelines-topo-150k.gdb")

CPU times: user 11.6 s, sys: 411 µs, total: 11.6 s
Wall time: 11.6 s


%%time
akl_roads = pd.clip(roads, AKL)

CPU times: user 23 s, sys: 0 ns, total: 23 s
Wall time: 23 s


# If a road has a highway number (hway_num not None), it's a highway/motorway
mway = akl_roads[~akl_roads.hway_num.isna()].copy()
mway


mway.name.value_counts().head(50).plot(kind="barh").invert_yaxis()


mway.hway_num.value_counts().head(50).plot(kind="barh").invert_yaxis()


ax = mway.to_crs(epsg=3857).plot()
ax.set_title("Auckland Region motorways")
ctx.add_basemap(ax)


%%time
mway_cube = make_geocube(
    vector_data=mway,
    measurements=["lane_count"],
    like=geocube, # Ensures dimensions match
    fill=0, # 0 works fine here, as every mway has at least one lane
)
mway_cube

CPU times: user 220 ms, sys: 0 ns, total: 220 ms
Wall time: 218 ms

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    lane_count   (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


mway_cube.lane_count.plot()
outfile = "output/mway.tif"
if not os.path.isfile(outfile):
    mway_cube.lane_count.rio.to_raster(outfile, dtype=np.byte)


src_ds = gdal.Open("output/mway.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/mway_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
mway_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(mway_dist.shape)
plt.imshow(mway_dist)
plt.title("Distance from motorways in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)

Text(0.5, 1.0, 'Distance (m)')


%%time
cbd = pop2018[pop2018.MB2020_V2_00 == "0433501"].copy()
cbd.geometry = cbd.geometry.buffer(1000)
cbd_cube = make_geocube(
    vector_data=cbd,
    like=geocube, # Ensures dimensions match
    fill=0
)
outfile = "output/cbd.tif"
cbd_cube.General_Electoral_Population.rio.to_raster(outfile, dtype=np.byte)
    
src_ds = gdal.Open("output/cbd.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/cbd_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
cbd_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(cbd_dist.shape)
plt.imshow(cbd_dist)
plt.title("Distance from CBD in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)
CPU times: user 298 ms, sys: 20.2 ms, total: 319 ms
Wall time: 315 ms

Text(0.5, 1.0, 'Distance (m)')


bounds = AKL.total_bounds.tolist()
bounds

[1703081.9789640256, 5870396.320936217, 1804839.668875325, 6002367.198185163]


zf = zipfile.ZipFile('input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip')
tiles = [file for file in zf.namelist() if file.endswith(".tif")]
tiles

['EJ.tif',
 'DM.tif',
 'EL.tif',
 'DL.tif',
 'DJ.tif',
 'FK.tif',
 'DK.tif',
 'EK.tif',
 'FL.tif']


tile_datasets = [rasterio.open(f'zip://input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip!{tile}') for tile in tiles]
DEM, transformation = rasterio.merge.merge(tile_datasets, bounds = bounds, res = (100,100), dtype=np.int16)


cbd_dist.shape

(1320, 1001)


# Ensure same shape
DEM = DEM[:,0:cbd_dist.shape[0],0:cbd_dist.shape[1]]
DEM.shape

(1, 1320, 1001)


print(np.nanmin(DEM), np.nanmean(DEM), np.nanmax(DEM), DEM.shape)
rasterio.plot.show(np.where(DEM>=0, DEM, np.nan), cmap='terrain')

-32767 -18180.83639996367 697 (1, 1320, 1001)

<AxesSubplot:>


meta = tile_datasets[0].meta
print(meta)
meta.update({
    "dtype": "int16",
    "height": DEM.shape[1],
    "width": DEM.shape[2],
    "transform": transformation
})
print(meta)
outfile = "output/slope.tif"
if not os.path.isfile(outfile):
    with rasterio.open(outfile, "w", **meta) as dest:
       dest.write(DEM)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -32767.0, 'width': 3028, 'height': 8192, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(8.0, 0.0, 1679712.0,
       0.0, -8.0, 5963776.0)}
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32767.0, 'width': 1001, 'height': 1320, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(100.0, 0.0, 1703081.9789640256,
       0.0, -100.0, 6002367.198185163)}


!ls -Ggh output

total 21M
-rw-r--r-- 1 1.3M Apr 15 18:26 cbd.tif
-rw-r--r-- 1 2.6M Apr 15 18:26 cbd_dist.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_1996.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2001.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2008.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2012.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2018.tif
-rw-r--r-- 1 1.3M Apr 15 14:37 mway.tif
-rw-r--r-- 1 2.6M Apr 15 18:26 mway_dist.tif
-rw-r--r-- 1 2.6M Apr 15 14:31 pop2013.tif
-rw-r--r-- 1 2.6M Apr 15 14:32 pop2018.tif
-rw-r--r-- 1 2.6M Apr 15 18:16 slope.tif

	Meshblock	MeshblockNumber	Population_Count_Usual_Resident_2013	Population_Count_Census_Night_2013	geometry
18549	MB 0712122	0712122	537	534	POLYGON ((1771009.117 5907877.430, 1770674.758...
10466	MB 0141001	0141001	228	222	POLYGON ((1749803.416 5970545.692, 1749808.725...
11552	MB 0181119	0181119	90	90	POLYGON ((1750734.441 5929622.619, 1750713.994...
13444	MB 0291101	0291101	159	159	POLYGON ((1748571.959 5912405.305, 1748576.611...
19198	MB 0746312	0746312	414	417	POLYGON ((1763875.903 5906059.287, 1763870.356...

	MB2020_V2_00	General_Electoral_Population	Maori_Electoral_Population	GED2020_V1_00	GED2020_V1_00_NAME	GED2020_V1_00_NAME_ASCII	MED2020_V1_00	MED2020_V1_00_NAME	MED2020_V1_00_NAME_ASCII	LAND_AREA_SQ_KM	AREA_SQ_KM	Shape_Length	geometry
16505	0560200	129	-999	025	Mt Roskill	Mt Roskill	3	Tāmaki Makaurau	Tamaki Makaurau	0.049186	0.049186	886.217648	POLYGON ((1755329.254 5914846.754, 1755307.113...
19082	0768007	183	27	048	Takanini	Takanini	3	Tāmaki Makaurau	Tamaki Makaurau	0.071362	0.071362	1286.891703	POLYGON ((1769995.098 5898877.583, 1770033.332...
12463	0181113	87	-999	055	Upper Harbour	Upper Harbour	5	Te Tai Tokerau	Te Tai Tokerau	0.019524	0.019524	679.854115	POLYGON ((1750649.804 5930127.249, 1750651.029...
15551	0460000	69	15	049	Tāmaki	Tamaki	3	Tāmaki Makaurau	Tamaki Makaurau	0.021341	0.021341	619.305796	POLYGON ((1762325.891 5918422.826, 1762328.052...
8155	4002769	117	18	039	Port Waikato	Port Waikato	1	Hauraki-Waikato	Hauraki-Waikato	4.632509	4.632509	14480.605616	POLYGON ((1751128.613 5882109.586, 1751136.090...

	t50_fid	name_ascii	macronated	name	hway_num	rna_sufi	lane_count	way_count	status	surface	geometry
512	100120610	KAIPARA COAST HIGHWAY	N	KAIPARA COAST HIGHWAY	16	3007739	2	None	None	sealed	LINESTRING (1732000.000 5944172.070, 1732048.5...
2933	3198057	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748581.508 5968975.145, 1748558.4...
2934	3198059	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748171.047 5971284.152, 1748129.9...
3320	3200754	PAERATA ROAD	N	PAERATA ROAD	22	3000260	2	None	None	sealed	LINESTRING (1767236.112 5888088.508, 1767244.3...
3324	3200792	UPPER HARBOUR MOTORWAY	N	UPPER HARBOUR MOTORWAY	18	3047073	4	None	None	sealed	LINESTRING (1747954.314 5927269.837, 1747970.0...
...	...	...	...	...	...	...	...	...	...	...	...
138240	100048291	AUCKLAND-WAIWERA MOTORWAY	N	AUCKLAND-WAIWERA MOTORWAY	1	3067966	7	None	None	sealed	LINESTRING (1755881.018 5922863.734, 1755886.4...
138301	100048432	AUCKLAND-HAMILTON MOTORWAY	N	AUCKLAND-HAMILTON MOTORWAY	1	3017109	1	None	None	sealed	LINESTRING (1765115.647 5909916.697, 1765092.7...
138337	100048532	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	4	None	None	sealed	LINESTRING (1748892.089 5949596.727, 1748892.0...
138369	100048589	PORT ALBERT ROAD	N	PORT ALBERT ROAD	16	3013274	2	None	None	sealed	LINESTRING (1734173.019 5980575.187, 1734175.6...
138680	100118365	SOUTH-WESTERN MOTORWAY	N	SOUTH-WESTERN MOTORWAY	20	3018532	4	None	None	sealed	LINESTRING (1760066.252 5908184.133, 1760043.9...

Land use¶

Population density¶

Distance from major roads¶

Distance from CBD¶

Slope¶

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
198093	High Producing Exotic Grassland	High Producing Exotic Grassland	High Producing Exotic Grassland	High Producing Exotic Grassland	High Producing Exotic Grassland	40	40	40	40	40	...	no	yes	yes	yes	yes	yes	Landcare Research	2014-06-30T00:00:00	lcdb1000405451	MULTIPOLYGON (((1722221.220 5681067.797, 17222...
14330	Deciduous Hardwoods	Deciduous Hardwoods	Deciduous Hardwoods	Deciduous Hardwoods	Deciduous Hardwoods	68	68	68	68	68	...	no	yes	yes	yes	yes	yes	Landcare Research	2004-06-30T00:00:00	lcdb2000174099	MULTIPOLYGON (((1322047.214 5038054.408, 13220...
5865	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	54	54	54	54	54	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000113289	MULTIPOLYGON (((1995230.346 5707153.821, 19952...
211672	Forest - Harvested	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	64	71	71	71	71	...	no	yes	yes	yes	yes	yes	Landcare Research	2019-12-01T00:00:00	lcdb2000505876	MULTIPOLYGON (((1573982.862 5259690.288, 15739...
132367	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	52	52	52	52	52	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000083894	MULTIPOLYGON (((2004675.083 5675750.092, 20046...

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
368793	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	69	69	69	69	69	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000165924	POLYGON ((1754495.703 5938405.984, 1754470.663...
176164	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	54	54	54	54	54	...	no	yes	yes	yes	yes	yes	Landcare Research	2004-06-30T00:00:00	lcdb1000122322	POLYGON ((1733140.709 5960181.994, 1733141.158...
306801	Built-up Area (settlement)	Built-up Area (settlement)	Built-up Area (settlement)	Built-up Area (settlement)	Built-up Area (settlement)	1	1	1	1	1	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000001559	POLYGON ((1773629.254 5899520.127, 1773671.247...
172842	High Producing Exotic Grassland	High Producing Exotic Grassland	Exotic Forest	Exotic Forest	Exotic Forest	40	40	71	71	71	...	no	yes	yes	yes	yes	yes	Landcare Research	2014-06-30T00:00:00	lcdb1000219851	POLYGON ((1714511.714 5957724.243, 1714527.942...
386662	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	69	69	69	69	69	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000164525	POLYGON ((1769735.325 5903460.892, 1769739.234...